import numpy
import scipy.sparse as sparse
import scipy.sparse.linalg as linalg
import matplotlib.pyplot as plt
import plotly.graph_objects as go
import pandas as pd
from ipywidgets import interact, SelectionSlider
import plotly.io as pio
import plotly.offline as py
pio.renderers.default = "notebook_connected"
py.offline.init_notebook_mode()
Numerical Analysis of 2D Diffusion-Reaction Equations Using an Implicit-Explicit Method: Solution and VisualizationΒΆ
Overview: This is a showcase of a hybrid implicit-explicit method to solve the diffusion-reaction PDE in 2-dimensions. It uses Plotly to provide a stroboscopic and interactive visualization of the solution.
Connection to data science: The PDE has indirect connections to time-series analysis and forecasting. The diffusion equation can be used to smooth time-series data, removing short-term fluctuations while preserving long-term trends. This is analogous to applying a low-pass filter. The diffusion helps in removing short-term fluctuations, and the reaction term can introduce additional effects such as trends or external influences.
Why use a hybrid method?: The diffusion term involves second spatial derivatives, which tend to amplify numerical errors in explicit methods like Forward Euler. The stability criterion for the diffusion term is $\Delta t=\frac{\Delta x^2}{4D}$, where $D$ is the diffusion coefficient [3]. This requires small time steps, making the method not time-efficient. By using an implicit method for the diffusion term and an explicit method for the reaction term, the hybrid approach is more time-efficient. β
The Implicit-Explicit Method:ΒΆ
The PDE is
\begin{aligned} u_t &= \sigma D_1 (\partial_{xx}+\partial_{xx})u + f(u, v) \\ v_t &= \sigma D_2 (\partial_{xx}+\partial_{xx})v + g(u, v) \end{aligned}
with the reaction terms
\begin{aligned} f(u,v) &= \alpha u (1 - \tau_1 v^2) + v (1 - \tau_2 u) \\ g(u,v) &= \beta v + \alpha \tau_1 u v^2 + u (\gamma + \tau_2 v). \end{aligned}
The following numerical method solves the PDE over the square domain $\{(x,y):-1\leq x,y\leq1\}$ with periodic boundary conditions.
Numerical method:
\begin{aligned} U^\ast &= U^n + \Delta t \sigma D_1 (\partial_{xx}+\partial_{xx}) U^\ast \\ V^\ast &= V^n + \Delta t \sigma D_2 (\partial_{xx}+\partial_{xx}) V^\ast \\ U^{n+1} &= U^\ast + \Delta t f(U^\ast, V^\ast) \\ V^{n+1} &= V^\ast + \Delta t g(U^\ast, V^\ast). \end{aligned}
The following function distretizes the spacial derivative using a five-point stencil in 2d and returns a sparse matrix reprsentation.
def laplacian_discretization(m):
"""Constructs a sparse matrix that discretizes the 2d Laplacian
Uses a five-point stencil and periodic boundary conditions.
"""
delta_x = 2.0 / (m + 1)
# Primary discretization
e = numpy.ones(m)
T = sparse.spdiags([e, -4.0 * e, e], [-1, 0, 1], m, m)
S = sparse.spdiags([e, e], [-1, 1], m, m)
I = sparse.eye(m)
A = sparse.kron(I, T) + sparse.kron(S, I)
# Construct periodic BCs
e = numpy.ones(m**2)
A_periodic = sparse.spdiags([e, e],[m - m**2, m**2 - m], m**2, m**2).tolil()
# Left & right BCs:
for i in range(m):
A_periodic[i * m, (i + 1) * m - 1] = 1.0
A_periodic[(i + 1) * m - 1, i * m] = 1.0
# Combine two matrices
A = A + A_periodic
A /= delta_x**2
A = A.todia()
return A
def f_reaction(U, V, sigma, tau_1, tau_2, alpha, beta, gamma):
return alpha * U * (1.0 - tau_1 * V**2) + V * (1.0 - tau_2 * U)
def g_reaction(U, V, sigma, tau_1, tau_2, alpha, beta, gamma):
return beta * V * (1.0 + alpha * tau_1 / beta * U * V) + U * (gamma + tau_2 * V)
def forward_euler_step(U, V, delta_t, A, sigma, f, g, D1=0.5, D2=1.0):
"""Take a single forward Euler step on the reaction-diffusion equation"""
U_new = U + delta_t * (sigma * D1 * A * U + f(U, V))
V_new = V + delta_t * (sigma * D2 * A * V + g(U, V))
return U_new, V_new
def backward_euler_diffusion_step(U, V, A, delta_t, sigma, D_1, D_2):
"""Take a single backward Euler step on the reaction-diffusion equation"""
U = linalg.spsolve((sparse.eye(A.shape[0]) - delta_t * sigma * D_1 * A), U)
V = linalg.spsolve((sparse.eye(A.shape[0]) - delta_t * sigma * D_2 * A), V)
return U, V
def imex_solver(sigma, tau_1, tau_2, alpha, beta, gamma, D_1, D_2):
# Alias reaction functions with the above parameters
f = lambda U, V: f_reaction(U, V, sigma, tau_1, tau_2, alpha, beta, gamma)
g = lambda U, V: g_reaction(U, V, sigma, tau_1, tau_2, alpha, beta, gamma)
# Set up grid
m = 150
delta_x = 2.0 / m
x = numpy.linspace(-1.0, 1.0, m)
y = numpy.linspace(-1.0, 1.0, m)
Y, X = numpy.meshgrid(y, x)
# Initial data
U = numpy.random.randn(m, m) / 2.0
V = numpy.random.randn(m, m) / 2.0
#fig = plt.figure()
#axes = fig.add_subplot(1, 1, 1, aspect='equal')
#plot = axes.pcolor(x, y, U, cmap=plt.get_cmap("viridis"))
#fig.colorbar(plot)
# Setup spatial discretization
U = U.reshape(-1)
V = V.reshape(-1)
A = laplacian_discretization(m)
# Time
t = 0.0
t_final = 30.0
delta_t = delta_x / (10.0 * sigma)
num_steps = int(numpy.round(t_final / delta_t))
data = []
data.append([x,y,U, t])
# Evolve in time
next_output_time = 0.0
for j in range(num_steps):
U, V = backward_euler_diffusion_step(U, V, A, delta_t, sigma, D_1, D_2)
U, V = forward_euler_step(U, V, delta_t, A, sigma, f, g)
t += delta_t
if t >= next_output_time:
next_output_time += 5.0
U_output = U.reshape((m, m))
#fig = plt.figure()
#axes = fig.add_subplot(1, 1, 1, aspect='equal')
#plot = axes.pcolor(x, y, U_output, cmap=plt.get_cmap("viridis"))
#fig.colorbar(plot)
#axes.set_title("t = %s" % t)
data.append([x,y,U_output, t])
#plt.show()
return data
# Parameters
data = imex_solver(sigma=0.0021, tau_1=3.5, tau_2=0, alpha=0.899, beta=-0.91, gamma=-0.899, D_1=0.5, D_2=1.0)
print('checkpoint')
checkpoint
# define time_intervals
time_intervals = []
for i in range(len(data)):
time_intervals.append(data[i][3])
def extract_data(time):
j = time_intervals.index(time)
z_data = data[j][2]
return z_data
def update_plot(time):
z_data = extract_data(time)
z_data = z_data.reshape((150, 150))
fig = go.Figure(data=[go.Surface(z=z_data)])
fig.update_traces(contours_z=dict(show=True, usecolormap=True,
highlightcolor="limegreen", project_z=True),)
fig.update_layout(title='Topographical 3-D Surface Plot of the Solution with Contours at t={}'.format(time), autosize=False,
scene_camera_eye=dict(x=1.7, y=1.3, z=0.64),
width=1000, height=700,
margin=dict(l=85, r=50, b=65, t=90),
)
fig.show()
for i in time_intervals:
update_plot(i)
# These widgets are for a Jupyter environment
#time_slider = SelectionSlider(options=time_intervals, value=time_intervals[5], description='Time')
#interact(update_plot, time=time_slider);
ReferencesΒΆ
[1] Mandli, K. T. (n.d.). Reaction-Diffusion Demo. In Numerical Methods for Partial Differential Equations [Jupyter Notebook]. GitHub. Retrieved July 27, 2024, from https://github.com/mandli/numerical-methods-pdes/blob/master/reaction-diffusion_demo.ipynb
[2] Ketcheson, D. (n.d.). Finite Difference Course. GitHub. Retrieved July 27, 2024, from https://github.com/ketch/finite-difference-course
[3] LeVeque, R. J. (2007). Finite Difference Methods for Ordinary and Partial Differential Equations: Steady State and Time Dependent Problems. Society for Industrial and Applied Mathematics (SIAM).